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Abstract. The thickness of the neutral hydrogen layer, 
coupled with the rotation curve, traces the outer dark 
matter potential. We estimate the amplitude of the flaring 
in spiral galaxies from a 3D model of the HI gas. Warps 
in particular are explicitly parametrized in the form of 
an harmonical density wave. Applying our method to the 
galaxy NGC 891, the only model that could fit the ob- 
servations, and in particular the HI at large height above 
the plane, includes a strong warp with a line of node al- 
most coinciding with the line of sight. This high-Z HI is 
not observed at the most extreme velocity channels, those 
corresponding to high rotational velocities. This is ac- 
counted for by the model, since orbits in the tilted planes 
are not circular, but elongated, with their minor axis in 
the galaxy plane. Their velocity on the major axis (i.e. at 
their maximal height above the plane) is then 30% less 
than in the plane. We finally connect the modelled verti- 
cal outer gaseous distribution to the dark matter through 
hydrodynamical and gravitational equations. Under the 
assumption of isotropy of the gaseous velocity dispersion, 
we conclude on a very flattened halo geometry for the 
galaxy NGC 891 (q w 0.2), while a vertical velocity dis- 
persion smaller that the radial one would lead to a less 
flattened Dark Matter Halo (q ps 0.4 - 0.5). Both results 
however suggests that dark matter is dissipative or has 
been strongly influenced by the gas dynamics. 

Key words: Galaxies: general - Galaxies: halos - Galax- 
ies: individual: NGC891 - Galaxies: ISM - Galaxies: kine- 
matics and dynamics - Cosmology: dark matter 



1. Introduction 

Knowing the 3D shape and extent of the Dark Matter 
Halos embedding galaxies (hereafter DMHs) is of partic- 
ular importance: it can constrain the nature of the dark 
matter itself, and in particular its dissipative character. 
Even in the frame of a given dark matter model, such 
as CDM, the shape of DMHs can constrain the galaxy 
formation scenarios, including gas infall (e.g. Dubinski & 



Carlberg 1991, Dubinski 1994). Little is known about the 
shape of DMHs, and in particular their flattening and ex- 
tent. Most of the best evidence of the existence of DMHs 
has been gained through HI rotation curves (see the re- 
view by Freeman 1993), and are limited to the HI disk. 
The dynamics of satellites (Zaritsky & White 1994) can 
extend further but with strong biases. Techniques based 
on gravitational lensing, virial masses or X-rays can reach 
constraints at much larger scales. In a recent paper, David 
et al (1995) observed a decrease of the mass-to-light ra- 
tio [M to t /Mi um ] between galaxies and clusters, suggesting 
that the dark matter in the universe is essentially confined 
around individual galaxies. 

Since the best tool to test DMHs, HI rotation curves, 
only provide azimutal constraints, we have to find a verti- 
cal signature of these halos on the plane of spiral galaxies 
(for which we have the most accurate rotation curves). The 
knowledge of both the HI velocity dispersion and plane 
thickness is a mean to test the shape of DMHs. Since it 
has been measured in face-on galaxies a constant HI ve- 
locity dispersion with radius (e.g. Dickey et al 1990), the 
idea is to use the flaring of the neutral hydrogen (i.e. radial 
increase of the thickness of the HI layer towards external 
parts), as a vertical tracer of the dark matter potential. 
The HI layer in an isolated exponential stellar disk is ex- 
pected to flare exponentially with radius; but the ampli- 
tude of the flaring is expected to be modified with the pres- 
ence of surrounding dark mass, depending on its geometry. 
The first serious attempt using the flaring of the gas lay- 
ers has been done by R. Oiling (1995) who concluded on 
a very pronounced flattening of the halo of NGC 4244 of 
(| =)q = 0.2t°;f (i.e. E5-E9), the large error being mostly 
due to the uncertainties on the velocity dispersion of the 
gas. Other techniques exist: the first one is to constrain 
the flattening using Polar Rings (i.e. two "orthogonal" ro- 
tation curves). Sackett and Sparke (1990) found using this 
way a constraint on NGC4650a 's DMH of E6-E7 , while 
more recently the same study placed a E5 constraint on 
A0136-0801's DMH (Sackett and Pogge, 1995). However, 
the method using Polar Rings is intrinsically uncertain, 
and the Polar Ring itself can totally change the DMH 
geometry, as shown by Combes & Arnaboldi (1995). A 
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related method consists in using the dynamics of a pre- 
cessing dusty disk: Steiman-Cameron et al (1992) inferred 
this way for the SO galaxy NGC4753 a quasi-spherical El 
DMH. Deriving the DMH geometry by the flaring law has 
however the advantage to be applicable to a larger sample 
of galaxies : every spiral galaxy with an inclination > 60° 
can be an acceptable candidate (Oiling 1995). 
The large range of values found until now for DMH Batten- 
ings (from El to E9) can be interpreted either in the way 
that the local universe is very heterogeneous, or, and this 
is more probable, that methods quoted above are domi- 
nated by errors and uncertainties. 

Determining the flaring of the gas in disk galaxies is a 
difficult task, due to projection uncertainties (galaxies are 
never perfectly edge-on) and due to the presence of warps, 
that deform the gas plane and thus modifies our percep- 
tion of its vertical distribution. We remark that in the lit- 
erature concerning the determination of the flaring of the 
gas in spiral galaxies, warps are never taken into account, 
so that the validity of the flaring curves obtained remains 
correct only for the inner unwarped disk, precisely where 
the dark matter halo is undetectable: a statistical analysis 
of the morphology of 167 rotation curves of spiral galaxies 
indeed shows that the dynamics inside the optical regions 
is not usually controlled by the dark component (Corradi 
& Capaccioli, 1990, and Freeman, 1993). A predictive and 
consistent theory of warps still does not exist (Binney, 
1992) even if we may have new hints about it (Masset & 
Tagger, 1996). We still do not have a general model of 
warps, an equation connecting the shape of the warp to 
the fundamental parameters of the galaxy, so that most 
of the time, the warp is ignored in models. The flaring 
curves given in the litterature are therefore overestimated 
(nearly all galaxies being more or less warped). The DMH 
flattening is consequently generally underestimated (see 
figure 10 of Olling's Thesis, 1995 and our Fig 15): halos 
are certainly even more "disk-like" if we model and sub- 
stract a warp. More generally and as far as we know, no 
attempt has been made to search different "geometries" 
than ellipsoidal. Combes and Arnaboldi (1995) suggested 
that dark matter could be ring-like, in the same category 
of geometry than the Pfcnnigcr & Combes (1994) model of 
dark matter which is a disk. A question arises finally: can 
we use "warps" as a supplement vertical signature ? The 
idea that one may take the shape of warps as another con- 
straint on the DMH geometry (Sparke & Casertano 1988) 
is more and more controversial. Tubbs & Sanders (1979) 
and Sparke & Casertano (1988) developped the idea that 
the shape of warps can help to constrain the DMH geom- 
etry. But this is extremely model-dependent. A spherical 
potential can maintain a warp forever (Tubbs & Sanders, 
1979) while a self-gravitating warp also solves the prob- 
lem of differential precession (Sparke & Casertano 1988, 
Pfenniger et al 1994). The HI thickness might thus be the 
more reliable vertical signature of the halo on the outer 
disk. 



2. Modeling the 3D distribution of the HI gas in 
spiral galaxies 

In this paper we develop a theoretical method to derive 
constraints on the DMH flattening from the observed HI 
thickness and we apply it to the nearby galaxy NGC 891, 
which is almost edge-on, and on which an abundant lit- 
erature exists. The HI data of this galaxy used here are 
from Rupen (1991). NGC 891 is classified as an Sb galaxy, 
member of the NGC 1023 group, at a distance of 10.0 Mpc 
(Rupen 1991, for H = 75 km/s/Mpc) and of systemic 
velocity 535 km/s. This galaxy is not obviously warped, 
and it has a gaseous large-scale asymmetry (southern ex- 
tension up to Rhi, south = 26.5 kpc), while in the north 
the gas stops surprisingly much before the optical radius 
(Rhi, north = 18.375 kpc while R opt = 22.3 kpc). Our sub- 
ject of interest is to determine if there is indeed a warp, 
then estimate the true HI thickness and amplitude of the 
flaring. The entire modeling is geometrical and kinemati- 
cal, so that the question of assuming any time-dependency 
of warps or the spiral arms do not raise here. We give our- 
selves a 3D HI model and adapt the parameters to the 
observations, we do not simulate an evolution of the sys- 
tem as does the N-Body technique. 

2.1. Modeling the HI infinitely thin disk from the inte- 
grated flux along the line of sight 

The gaseous surface density can be represented by the real 
part of the Fourier development: 

oc 

S fl (i2,0) - A (R) + J2^n(R)e l{k - R - m{e - e - )) 
1=1 

For quasi edge-on galaxies, because of the projection, the 
O information vanishes, and spiral arms, in the best case, 
can be rebuilt with spectral cube modeling (section 4), but 
not found with the technique used here. The technique 
is to consider the zero moment map of a nearly edge-on 
galaxy (i.e. the integrated HI flux). The r.m.s. noise in the 
channel maps of NGC 891, associated with the synthesized 
beam of FWHM = 20" x 20", is 0.780 mJy/beam. The 
velocity resolution, after Hanning smoothing, is 20.7 km/s. 
The integrated HI flux is obtained in making a slice ex- 
actly along the major axis of the galaxy providing a curve: 
g(X) (full line in Fig 1 and Fig 2) , where X is the radius 
along the major axis (corresponding to the X coordinate 
of the natural cartesian frame of the galaxy). 
This curve, g(X), is connected to the radial gaseous distri- 
bution E 9 (R) through an integral that can be inverted as- 
suming low optical depth (Abel inversion formula), yield- 
ing: 

E (R) - 1 f°°MX) dX 

The results are shown (dashed curve) in Fig 1 and 
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Fig. 1. The full line represents the observational integrated HI 
flux along the northern part of NGC 891 from r.m.s corrected 
zero moment map (in Jy/beam), while the dashed line is the 
surface gaseous density inverted from the integrated HI flux 
assuming low optical depth (in Mq/pc 2 ). The dot-dashed-dot 
gaussian curve represents the best fit of the axisymmetric 
gaseous density (m — 0), taking the southern part into ac- 
count. 



Fig 2. We note a good symmetry of the two curves inside 
Riii,north so that at first approximation we may write 
E S (J?, 9) = T, g (R) = A {R). We fitted the gaseous ra- 
dial distribution quite well with a gaussian law of charac- 
teristics given below. There are two incertainties in this 
determination. First, the density formula (1) give nega- 
tive gaseous densities in the very central parts (R < 2 
kpc), since the derivative of the integrated flux curve is 
positive, such a behaviour of the g(X) curve being due 
to an HI depletion in the center of the NGC 891 galaxy, 
which is quite often seen in the central parts of spiral 
galaxies. Second, the integrated flux has an intrinsic error 
due the non-infinite resolution, and we could observe that 
the density inversion formula was highly sensitive to lo- 
cal variations of the integrated flux curve. Assuming that 
the HI is distributed in a disk of radius Rhi,o consists 
in replacing in the right term of (1) oo by Rhi,o- It was 
interesting to find that the resulting radial density of the 
HI gas remains practically the same with or without this 
assumption (for Rhi,o = 19 kpc). This suggests that the 



Fig. 2. Same as Fig 1 but for the southern part of NGC 891. 
We can observe a relatively good symmetry with the northern 
part until 5.5 arcmin (fitted with the gaussian curve), confirm- 
ing the predominance of the axisymmetric mode. 



high-Z gas surrounding the galaxy, which is probably lo- 
cated at large radii, is no more confined in the thin central 
HI disk, finally suggesting the presence of a warp along the 
line of sight, which is also the direction of integration of 
the HI flux. The gaussian law fitting axisymmetrically the 
northern and southern part of NGC 891 is: 



MR) = 



2.354p 0;i// 



Qg 



2n 



exp(- 



(R 



2g 



(2) 



with R g = 9.5 kpc, g g = 3.8 kpc and po.Hi = 9.6 Mq/pc 2 . 
This fit of the axisymmetric radial gaseous distribution is 
quite sufficient for the rotation curve fit (next section), but 
for spectral cube modeling, we kept integrally the dashed 
curve representing the "real" Y, g (R), through stored data. 
We finally calculate the relation between the total HI mass 
and a gaussian surface density: 



M HI W 2TTp ,HlQgRg(^+2^ 



This gives for M HI = 4.2 10 9 M e /pc 2 (Rupen 1991 for 
h = 0.75) the value po,Hi = 8.0 Mq/pc 2 , in good agree- 
ment with our value of Po,hi = 9.6 Mq/pc 2 . 
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2. 2. Modeling the warping of the disk 

The warp is modelled as a vertical density wave of equa- 
tion (m=l): 



Z warp (R,e) = T(R)cos(Q- e w ) 



(3) 



The amplitude of this vertical wave, T(R), is linearly 
parametrized but we also tested a law in R 2 (paraboloid): 



T{R < R w ) 
T(R > R w ) 
T(R > R w ) 





an 



a 



par 



(R- 
(R 2 



Rw) 



This parametrization corresponds to a flat disk's plane 
until a radius R w , where the warp starts, and after fol- 
lows the shape of the T function. We primarily included a 
global wave number for the warp wave (a k w R term in the 
cosine) but we finally neglected it for two reasons: firstly 
because the line of nodes of the warps are generally ob- 
served to be almost straight until large radii (Briggs, 1990) 
yielding that k w = 0, and secondly because a twist of the 
line of nodes (k w ^ 0), starting after some radius, is any- 
way hard to detect with spectral cube modeling (section 
4) . Simulating modes of warps embedded in dark matter 
halos, Sparke & Casertano (1988) frequently found linear 
solutions for the shape of the T functions. Briggs (1990) 
determined the tip angles (i.e. T(R)) of warps of a sample 
of 12 galaxies. Nearly all the features derived follow glob- 
ally a linear or paraboloidal shape. Warps are generally 
observed to raise up at 3-5 exponential scale lengths of 
the light distribution, which corresponds more or less to 
the location of the optical radius. We show in Fig 11a 3D 
model of NGC 891 which revealed to possess a warp of 
high amplitude, beginning at 2.5 exponential scale length 
(12 kpc). 

2. 3. Modeling the HI vertical distribution 

It is expected from hydrostatic equilibrium that the HI 
layer is vertically distributed under a gaussian law only if 
the gas layer is much thinner than the stellar scale height. 
For the case of NGC 891, the exponential stellar disk has 
a vertical scale height hz ~ 1-0 kpc, and we determine 
in a further section, using our 3D spectral cube modeling, 
that the height of the inner gaseous layer H is about 0.1 
kpc so that we can assume a gaussian vertical distribution 
for the HI gas: 



p gas oc exp(- 



Z 2 



2Z (R) 2 



(4) 



If we write like Oiling (1995) the Hydrostatic Equilibrium 
of the HI gas, under the assumption of Z-isothcrmality 
and axisymmetry, we obtain: 



2 d\np HI 



galaxy 

dZ 



INNER FLARING OF NGC 891 




R (kpc) 



Fig. 3. NGC 891 flaring for a gaussian vertical distribution 
from the equation of hydrostatic equilibrium, using the rota- 
tion curve fit section results. This curve is invariant for Z be- 
tween 0.1 and 0.5 kpc. Rwarp is determined in the spectral 
cube modeling section. 



so that inserting (4) in (5) yields: 



Z (R) 



Z 



9*„ 



dZ 



-(R, Z) 



(6) 



Equation (6) is correct for small Z, and R smaller than 
R Wl since after this radius, the gas rises up in the warp. 
Equation (4) is then modified. Its form is also dependent 
on the dark matter model adopted: a disk-like dark mat- 
ter distribution, following the HI layer, settles the gas in 
a potential well and the resolution of equation of the hy- 
drostatic equilibrium will yield a solution in sech-2, while 
a spherical dark matter halo will provide a vertical force 
proportionnal to Z and thus a gaussian solution for the 
equation. But a sech-2 function and a gaussian have com- 
parable behaviours for Z smaller than 1 kpc, so that in 
practice they cannot be separated. We thus adopted a 
gaussian representation: 



(5) Pgas oc exp(- 



(Z — Z warp (R, < d)) z 
2Z (R) 2 



(7) 



J.-F. Becquaert and F. Combes: The 3D Geometry of Dark Matter Halos 



5 



We derive the galaxy potential of (6) from the rotation 
curve fit of next section. We determine using Fig 3: 

Z (R < 12 kpc) = 0.08 + 0.014R [kpc] 



o 
o 

CO 



NGC 891 Rotation Curve Fit 



We determine further (spectral cube modeling section) 
that this galaxy is warped relatively early in radius (be- 
ginning at about 12 kpc from the center of the galaxy, sim- 
ilarly to the Milky Way) so that an outer flaring cannot 
be determined only with a rotation curve fit. We observe 
that the inner flaring obtained above is almost linear as 
expected: many authors, (Pfcnnigcr et al, 1994 or Oiling, 
1995) predicted such a behaviour. This leads us to consider 
for the spectral cube modeling section a warped gaussian 
HI vertical distribution with a flaring of the form: 



Z {R < R f ) 
Z (R > R f ) 



H 
H 



CiR 



UR-Rf) 



(8) 



where (i is the inner flaring coefficient and ( 2 the outer 
flaring coefficient. The £ coefficients are the slope of the 
flaring, they are dimensionless. The (2 coefficient is nat- 
urally the most sensitive to the DMH geometry. This 
parameter can not be satisfactorily approached without 
modeling in 3 dimensions the warped neutral hydrogen 
layer. 

3. Fit of the rotation curve 

NGC 891 having a quasi edge-on orientation, and because 
of HI deficiency in its center, the HI rotation curve is in- 
nacurate inside 8 kpc. The HI points are from Sancisi & 
Allen (1979), the relative error being of about 5 % af- 
ter 8 kpc, in the flat region. A multi-component fit, al- 
though the method gives non-unique solutions, requires 
at least sufficiently good data, so that we considered for 
radii smaller than 8 kpc the CO curve provided by Garcia- 
Burillo et al (1992). NGC891 is expected to have a fast ro- 
tating nuclear disk (quite similarly to the Milky Way, see 
Dame et al, 1987), explaining why the CO rotation curve 
steepens to 250 km/s in 500 pc to fall at 225 km/s and stay 
constant up to 5 % at 225 km/s. We consider hereafter five 
components: a Toomrc-Kuzmin nuclear disk, a Plummer 
bulge, a thick stellar disk, a gaseous contribution calcu- 
lated for the previously derived radial HI surface density, 
and an ellipsoidal dark matter halo. The nuclear Toomre 
disk has a mass M nuc i ear and a scale- length s nuc i ear . The 
Bulge spherical potential is classically 



*b(R, Z) 



-GM b 



^JR? + Z 2 + s 



and its rotational contribution is straightforward. The 
stellar disk is a thick exponential disk of density (van der 
Kruit, 1981): 

R Z 
p±(R<R op t,Z) = p ,*cxp(- — ) cxp(- — ) 

Hp 



o 
o 

cv 



o 
> 



o 
o 



Stellar disk - 



Nuclear disk_. 



10 15 
R (kpc) 



20 



25 



Fig. 4. Multicomponent fit of the rotation curve. Open squares 
are CO data from Garcia-Burillo et al (1992) while little full 
squares are HI data from Sancisi & Allen (1979). For this fit 
the dark component has a flattening of q = 0.5 and a core 
radius of 14 kpc. The HI contribution, negligible, is not drawn. 



p*(R>R opU Z) = 0.0 

The value of hz and Iir are kept fixed since several previ- 
ous studies have settled them with confidence (see Table 
1). The disk potential is then (Sackett & Sparke 1990): 



f 

Jo 



Jo{yR) [yZ d e 



-\z\/z d 



-v\z\] 



(1 + ^)3/2(^-1) 



-dy 



We derive the contribution of the stellar disk to the rota- 
tion curve: 



-GM+R 



f 

Jo 



vMvR) 



(1 + 42,2)3/2(^+1) 



dy 



hi 



We simulated also the equations of Casertano (1983) who 
includes explicitly the truncation of the stellar disk. But 
considering that the HI data stops before the optical ra- 
dius in the north, and that some HI in the south is prob- 
ably missing along the line of sight, truncation effects on 
the rotation curve are not expected to be visible, so that 
we adopted the formula of de Zeeuw & Pfcnnigcr (refound 
by Sackett & Sparke, 1990). 
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Fig. 5. The impact of modifying the dark matter density model 
on the rotation curve fit. Values of the parameters of the dark 
matter halo component are kept fixed, while only the power p 
changes. 



We finally consider a dark matter component of isother- 
mal, pseudo-ellipsoidal density: 



Pdmh 



Pq.dmh 



+ W, + W^) 



This distribution generates the potential: 
Ydmh(R, Z) = 2irGpo.DMHqRc 



f 

Jo 



Hn[l + fk^ + ^ 2 )] 



e 2 x 2 + 1 



dx 



where e = \J\ — q 2 . The halo contribution to the rotation 
curve is therefore: 



V DMH = ^Gpo,DMHqR 2 Rl 



f 

Jo 



x 2 dx 



l R 2 (e 2 x 2 + l) 2 + (e 2 x 2 + l)x 2 R 2 
which gives an asymptotic value, at large radii: 
2 „ n p2 garccosq 

v DMH,oo — 47TCrp , DMHn. c —f=== 

v 1 - q 



Fig. 6. The central density of the DMH as a function of its 
flattening, as found from the rotation curve fit. 



An interesting particular case of dark matter density 
model is: 



Pdmh — 



Pa, DMH 



z 1 



(1 + ^2 + R 2 tl 



(P>0) 



The corresponding formula for the potential and rotation 
curve contributions are given in appendix C. The intro- 
duction of the power parameter p (we retrieve the normal 
dark matter density model for p = 1) is justified by the 
work of Lake & Feinswog (1989), who considered the pos- 
sibility that Pdmh(R, Z) is not necessarily in R~ 2 as it 
is currently assumed, the flatness of the rotation curves 
being also achievable for profiles in R~ 3 and even R~ A 
(corresponding respectively to p — | and p — 2). We sim- 
ulated both values of p at the end of this section, in order 
to estimate the implication on the amount of dark matter 
required (Fig 5). 

Taking the gaussian fit of HI radial distribution, we derive 
using formula 2-160 of BT87 the HI contribution: 



J HI 



f-OO 

-Gp 0t HiR / xJi{xR) 
Jo 



exp(- 



{R-R g f 

2Ql 



)Jo{xy)ydy]dx 



J.-F. Becquaert and F. Combes: The 3D Geometry of Dark Matter Halos 



7 



In fact, the HI contribution to the rotation curve revealed 
to be negligible: in the outer parts, it is of the same order 
than the nuclear disk contribution. Inserting the HI con- 
tribution modifies the dark mass required of about 2 %, 
less than the error made in ignoring the stellar truncation. 
The halo parameter R c can be somehow constrained us- 
ing the [M/L] ratio of the stellar disk. Several simulations 
were made for several halo core radii: in each case the flat- 
ness of the rotation curve could be obtained, but for R c 
values lower than 12 kpc, the disk mass had to be more 
and more severely decreased. The disk mass required for 
R c = 14 kpc is M* = 7.7 1O 1O M yielding a mass-to-light 
ratio in blue color of [M/L]* : b — 7.1 M Q /L Q , in excellent 
agreement with the value of 7.0 found by van der Kruit 
(1981). For R c = 8 kpc, keeping the flatness of the curve 
requires the disk mass to be M* = 5.4 1O 1O M providing 
a mass-to-light ratio of 4.9 which is far too low (Bottema 
et al 1991). So that we settle the halo core radius at about 
14 ± 2 kpc. 

The mass of dark matter enclosed in an oblate ellipsoid of 
half major axis d and flattening q is given by: 

M H (<d) = 4:irp 0tDM HRlq[d- i? c arctan(-J-)] 

Rc 

Considering the dark mass enclosed within 30 kpc, in a 
semi-flattened ellipsoid of q=0.5, of core radius R c = 14 
kpc, we calculate that M H = 17 1O 1O M . 
The simulation of the generalized dark matter densities 
with p = | and p = 2 modifies drastically the results. 
We show in Fig 5 the alteration of the dark matter con- 
tribution for the same halo parameters, but for different 
values of p. Non-isothermal values of p (i.e. p > 1) require 
much greater dark matter central densities, and modifies 
the stellar disk mass severely. We keep p = 1 for the rest 
of the article. We finally remark that fitting the HI data 
only (i.e. without considering the inner CO data), leaves 
unchanged the dark matter component parameters. The 
results of our best fit (Fig 4) are quoted in Table 1. We 
plot in Fig 6 the relation between po,dmh and q. 

4. Numerical method 

4--1- The code 

The construction of synthetic spectral cubes requires the 
use of many parameters, among which six are considered 
as free parameters: the warp parameters a w ,R w and the 
flaring parameters Rf, Hq, (1,(2- 

A rapid code is thus necessary. We separate the density 
model building from the kinematics (velocity sampling). 
To build the density model, and to project it, we choose 
to distribute randomly particles in space: this technique 
is much more efficient (in terms of CPU time) than a 3D 
cube construction, for a given spatial resolution; this is 
due to the fact that the 3D cube has a very little volume 
filling factor, therefore the number of particules required 



Table 1. Results of the rotation curve fit. 



quantity fitting value 





Snuclear 






0.4 kpc 




M nuc lear 




— 


1.5 10 Mq 




SB 




= 


2.5 kpc 




M B 






3.2 lO lo M 




h z 




= 


0.99 kpc 




h R 




— 


4.9 kpc 








— 


7.7 10 Mq 




[M/L]*,b 




— 


7.1 Mq/Lq 




R c 






14 kpc 
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are much lower than the number of cube cells, and various 
3D rotations are then much quicker. 

4.1.1. Rendering the HI 3D geometry 

The computational form of the model density is: 

PHI (R,e,z) = x g (R) cxp(- (z 'f—^ e))2 ) 

where T, g (R) has been derived from section 2.1, approxi- 
mated by the equation (2), Z warp (R, <d) is the warp wave 
described by equation (3) and the flaring of the gas is in- 
cluded in Zo(R) through the linear law of equation (8). 
We generate N particles [x(i),y(i), z(i);i € 1,...,N] from 
the density model above, on which we operate a spatial 
rotation $ 3D = $1(1, j) o ^(PA, k), where "I" is the incli- 
nation of the galaxy and "P. A." its position angle, thus 
projecting the set of particles on the sky plane with the 
orientation angles of the real galaxy. 

4.1.2. Spectral sampling and synthetic spectral cube 

A velocity is associated to every particle in the cube. The 
rotational velocity V rot (R) is interpolated from the previ- 
ously fitted rotation curve. We project the velocities onto 
the plane of the sky to get radial velocities: 

Vios = V sys + V rot sin I COS 9 

Having a set of particles and the associated radial veloc- 
ities, we now sample this 4 dimensions space [x, y, z, V] 
according to our original NGC 891 spectral cube. This 
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sampling is done according to a central (systemic) veloc- 
ity and a fixed increment AV = \J AV^ han + a HI to get 
channel maps (this assumes an isotropic dispersion a hi). 
For our test-galaxy AV c h an is 20.7 km/s. The procedure 
of the construction is the following. For the Nst channel 
map Vjv = V sys + N x AV c han, we fill a temporary pixel 
cube with the spatially rotated particles in multiplying 
this geometric cube by a gaussian in Vjv, centered on our 
previously calculated Vi os , and of dispersion AV". We fi- 
nally sums the flux perpendicularly to the plane of the 
sky to obtain a channel map (the Nst). We finally loop 
over N (in our case from 5 to 25, 15 corresponding the 
systemic velocity) to get an entire modelled spectral cube 
[X, Y, Vios] - The last operation consists in convolving every 
pixel of the spectral cube by gaussians of FWHM corre- 
sponding to the synthetized beam of the VLA observations 
(20" x 20" for NGC 891). 

4.1.3. Application of the code 

The random launching of particles in space takes negligi- 
ble CPU, in front of the kinematics part. The dimensions 
of the cube is 300x300xN c ^ an , yielding a spatial resolution 
of 1 pixels 3 arcseconds (=147 parsecs), but the technique 
allows to build channels one by one, and only some of the 
total number of channels (N c ^ a „=30) are computed for 
the whole set of free parameters. In practice, 6 channel 
maps (whatever the number of particles) are sampled and 
beam-smoothed from one spectral cube in about 30 sec- 
onds on one processor of a CONVEX 3440. 
Spatial resolution is the dominant parameter to the con- 
trol of CPU time. The observed HI cube of NGC 891 
displays 30 channels, 24 of them having significant flux 
(channels 2 to 25), the center being the systemic velocity 
channel 15 (see Rupen 1991 for the mosaic figure of the 30 
channels). To restrict at maximum the CPU, we modelled 
every two channels. The inclination of this galaxy has been 
estimated several times in the literature, leading to a firm 
lower limit / > 88.6° (Rupen, 1991). The inclination is 
thus kept constant at 89.0 degrees since the uncertainty 
on its value is too small (±0.3 degree) to assign I to be 
a free parameter. In the same way, the position angle was 
fixed to 23.5 degrees (Garcia-Burillo et al 1992). Van der 
Kruit (1981), using Velocity-Position diagrams of models, 
concluded to a constant velocity dispersion of the gas (as- 
sumed isotropic) of 10.0 km/s. 

4.1.4. The experiments 

A library of models exploring the free parameters are com- 
pared to the HI observations (in figures 7, 8 and 9). 
With linear contours centered on the half-intensity of the 
channels, the spectral cube of NGC 891 (whether south or 
north) do not show any feature resembling to what Oiling 
nicely calls " butterfly wings" , which are characteristic of 
HI flaring and warp. To make weak structures at large 



height above the plane emerge, contours have to be ad- 
equately selected, as in Fig 7, which presents flaring and 
warped HI layers. 

Our first idea was that the surrounding high-Z gas (the 
low intensity features surrounding the disk up to several 
kpc height, visible in the channel maps "north NGC 891" 
of Fig 7, Fig 8 and Fig 9) could be interpreted as a warp 
in the line of sight, and so we fixed the value 8™ = 90°. 
We used to test this configuration the northern part of the 
spectral cube since no HI extension along the major axis 
is required. 

The radial distribution of the HI gas is assumed to be the 
dashed curve found in section 2.1. The gaussian model 
approximating this dashed curve is used later, for the de- 
tection of a spiral pattern; from that the axisymmetric 
mode is globally determined. We see in Fig 7, Fig 8, Fig 9 
samples of our intensive construction of spectral cubes. 
We vary nearly all possible parameters, intending to show 
what transformations on the 2D flux distribution corre- 
spond to the variation of each parameter. The first seven 
synthetic cubes shown in Fig 7 are presented to observe 
the 2D flux transformations that take place for extreme 
values of parameters. We do not search for this first se- 
quence to obtain cubes that are close to observations, but 
rather to build a library of models. 

xl: a Hn = 0, Ci = 0, ( 2 =0,H = 0.05 kpc. 

x2: a lm = 0, Ci = 0, C2 = 0, H = 0.9 kpc. This x2 
synthetic cube displays a too high HI constant thickness. 
By comparison with the observed cube, the average value 
is more towards low values (H ~ 0.1 — 0.5 kpc). 

x3: a lin (R w = 5 kpc) = 0.4, Ci = 0, C2 = 0, H = 0.4 
kpc. This synthetic cube displays a warp of too high am- 
plitude, generated too early. We clearly see that the HI 
disk flux of the observations is totally different. 

x4: a Un (R w = 1 kpc) = 0.2, Ci = 0, C2 = 0. The effect of 
decreasing the amplitude of the warp, makes its "wings" 
to join. The double "symmetric" holes of the channel 15 
are sufficient to rule out this model. 

x5: au n (R w = 3 kpc) = 0.08, (1 = 0, (2 = 0. For very 
low amplitudes of the warp, we reconstruct the observed 
HI disk: in fact, there is no more high-Z gas surrounding 
the disk, as is observed. 

x6: au n (R w — 10 kpc) — 1., (1 = 0, C2 = 0. We show 
here a warp generated at 10 kpc, further than before. We 
observe the emergence of high-Z gas surrounding the cen- 
tral HI disk. We also notice that R w is easy to locate in the 
way that the gas after R w is displaced, inducing a "trunca- 
ture" of the outer HI central disk. In particular, the 25th 
channel is very sensitive to that modification. We conclude 
from comparison with observations that the warp is prob- 




Fig. 7. Library of models: we present here 7 synthetic spectral cubes xl to x7 (to be compared with the original spectral cube 
"north NGC 891" displayed as the 8th cube). The same contours are used for every cube: 0.002:0.024:0.002 and 0.036 Jy/beam. 
The channels presented here are the systemic channel 15 (535 km/s), and the channels 17 (575 km/s), 19 (620 km/s), 21 (660 
km/s), 23 (705 km/s), 25 (745 km/s). 



ably generated a little bit further, around R — 12 kpc. 

x7: aii n = 0, £i = 0.15, (2 = 0, Hq = 0.1. We observe now 
the location of an exagerated inner flaring. This inner flar- 
ing of 0.15 is about ten times the inner flaring expected 
with the hydrostatic equilibrium of the gas (section 2.3). 

x8: a Hn = 0, Ci = 0, C2 = 0.6, R f = 12, H = 0.1. 
Testing the outer flaring coefficient at an abnormal value 
of 0.6. We observe that there are clear 2D effects induced. 

We now try to combine these results to build a sequence 
of cubes closer to the observed one. 

x9: a Un (R w = 10) = 1, Ci = 0.02, C2 = 0.02, H = 0.4. 
A flaring warp. We kept the warp parameters of cube x6, 
in adding a flaring of the gas. We notice that the flaring 
densifies the contours following the warp along the line of 
sight, so that even not very well constrained intrinsically, 
the flaring can be approached. 

The main problem when we observe the cube x9 in com- 



parison to the observed cube is that the extreme channel 
maps (23 and above all 25) are in every of our models sur- 
rounded by high-Z gas, while they are not in the observed 
cube. The observed channel map 25 is very thin, so that it 
appears necessary that the particles in the warp rotate at 
a different speed than the HI in the disk plane. This was 
already noticed by Swaters et al (1997), and interpreted 
as non-cylindrical rotation of the HI disk, due to energetic 
star-formation. To estimate quantitativally the amount of 
non-cylindrical rotation, we adopt a simple kinematical 
model: we assign a rotational velocity of high-Z particules 
lower than in the plane by a free factor (between 1 and 
2), for heights greater than the inner disk HI thickness, 
typically 150 parsecs. 

xlO: Simulating high-Z particles as rotating at Vq/2, where 
Vo is the "plateau" speed of the rotation curve in the plane 
(i.e. 225 km/s). The parameters for this cube are: 
a Un {R w = 9) = 1.3, Ci - 0.02, C2 = 0.02, H = 0.4. 
We observe that the high-Z gas concentrates around the 
center channels 15, 17, 19 leaving the extreme channel 23, 
25 without high-Z gas, which is one of the characteristics 




Fig. 8. We present 7 synthetic spectral cubes x8 to xl4 "converging" towards our best fits xl3 and xl4 (to be compared with 
the corresponding channels "north NGC 891"). The same contours are used in every cube: 0.002:0.024:0.002 0.036 Jy/beam. 



of the observed spectral cube. 

xll: same parameters but now with a parabolic warp in- 
stead of linear, of same maximum height as the previous 
cube xlO. We observe by comparison that a linear warp 
makes " round" waves instead of " square" waves obtained 
with a r function in R 2 . A linear warp seems closer to 
the observations. We also can appreciate that the warp is 
generated too early, truncating the end of the central HI 
disk (visible in particular for the 25th channel). 

xl2: same as xlO, but with V warp = Vb/1.3. We see that 
it is possible to make the 25th channel very thin, as in 
the observed cube. A velocity ratio of 1.3 shows up to be 
closer to the observations than the velocity ratio of 2 of 
the previous model cube since it allows the 23rd channel 
to have high-Z gas. 

xl3: A more accurate cube with V warp = V /1.25. We 
have here the definitive values of the geometric parame- 
ters: 

a Un {R w = 12) = 0.8, Ci = 0.02, C2 = 0.02, H = 0.1 
xl4: We refine the kinematical description of the warp 



behaviour detected. Instead of considering that the warp 
cylindrically rotates, we test now a vertical velocity gra- 
dient with an empirical " exponential" law. The model be- 
comes: 

Vios = V sys + V rot sin/cosecxp(-0.1|Z|) 

Regarding the hypothesis of the model (as the decoupling 
R-Z of the model) and the fact that the radial distribu- 
tion of the gas in section 2.1. is not perfect, we consider 
as successful (with the 2 last synthetic spectral cubes xl3 
and xl4) our goal to rebuild the observations. 

We try now to simulate the non-symmetric features ob- 
served in NGC 891 cube. We consider the sum of the R- 
gaussian axisymmetric mode Ao(R) determined in section 
2.1 and an m — 2 harmonic term. We generate the spi- 
ral after the HI central depletion of the galaxy, at about 
2 kpc, let it develop and make it decrease in l/R after a 
radius R sp i r ■ The formulation in the model is the following 

Y, g (R < 2 kpc) = 

H g {R < R spir ) = A {) {R) + A 2 cos{k 2 R-2{Q-Q 2 )) 
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Fig. 9. Spiral cubes: 6 synthetic spectral cubes representing an m = 2 spiral mode simulated alone, with a fixed wave number 
(k s = 0.05) and a varying line of nodes (from 0° to 150°). The height reached by the spiral has been a bit amplified in settling 
Rapir at 15 kpc, that is 3 kpc more than R w . We suppose that R sp ir is in reality of order R w . The seventh cube (x21) is a 
combination of a gaussian Ao(R) and an m = 2 harmonical mode with 02 = 30°. The same contours are used in every plot: 
0.002:0.024:0.002 0.036 Jy/beam. 



> Rspir) — Aq(R) H - ^-^ — — — 

where hi is the global wave number, O2 the line of nodes of 
the spiral pattern, and where A2 is of small magnitude in 
comparison to Aq(R). Computing alone (i.e. Aq(R) = 0) 
the 771 — 2 term in our spectral cube code (Fig 9), we 
could associate the external "torsion" of the gas, unex- 
plained with our symmetric models, to an m = 2 warped 
spiral. We kept the flaring warp parameters found in the 
axisymmetric mode construction. The velocity of the spi- 
ral arms is calculated in Appendix B and is implemented 
in the code. Synthetic cubes xl5 to x20 show an 777 = 2 
spiral with k 2 = 0.05 and regularly varying line of nodes 
(from 0° to 180°). It is important to notice that R sp i r is 
found greater than R w (of about 3 kpc), explaining the 
asymmetric features of the high-Z gas: the spiral is raised 
by the warp before its density vanishes, so that a non- 
negligible quantity of matter is asymmetrically placed in 
the line of sight, at heights corresponding to the warp 
wave. 



x21: This synthetic cube is the superposition of the ax- 
isymmetric gaussian density found in section 2.1., with an 
777 = 2 spiral (k 2 = 0.05, 62 = 30°). We observe unam- 
biguously that a spiral can reproduce the observed asym- 
metry. It is noteworthy that the gaussian HI density is only 
an approximation in front of the "real" £#/(i?) stored in 
buffers: the 25th channel, as an example, is less well fitted 
(see xl4 for comparison). 

4.1.5. Conclusion on the HI 3D distribution according to 
our model 

Several hundreds of synthetic cubes under the configura- 
tion of a warp along the line of sight (Q w = 90°) were 
finally built. Most parameters induce typical and particu- 
lar 2D effects: 

1- a flaring warp along the line of sight makes "butterfly- 
shapes" that can be associated to the observed high-Z 
features of the gas. 

2- a flaring also thickens regularly the central HI disk, a 
behaviour which is clear for the 23rd channel of NGC 891. 
Combinations step by step of the parameters succeeded 
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Fig. 11. The 3D HI isodensity contour of NGC 891 from a synthesis of our results. The radial distribution comes from 
section 2.1 while the vertical gaussian distribution has the parameters found in the spectral cube modeling: Ho — 0.1 kpc, 
Ci = C2 = C = 0.02. The warp appears clearly along the Y axis accordingly to 9„ = 90°. 



to reproduce most of the features observed in the origi- 
nal spectral cube of NGC 891. The radius where the warp 
begins, for instance, is very sensitive: we are able to set 
it with a good confidence at 12 ± 0.3 kpc. The amplitude 
of the warp, au n , is about 0.7 to 0.8, depending on the 
true pattern of the m = 2 spiral. The most important 
result is that we found possible to constrain the flaring 
(see synthetic cubes with varying (1 , £2 coefficients) . Let 
us summarize the other deductions: 

1- The flaring found after modeling in 3D is of the same 
order than the value estimated using the potential derived 
from the results of the rotation curve fit section (section 
2.3). 

2- Modeling with success the gas in 3D requires the use 
of every channel. The central channels were needed in our 
case to constrain the amplitude of the warp, and the ex- 
treme channels 23 and 25 were needed to detect the kine- 
matical behaviour of the warp wave. In particular, the 
channels 15 and 23 showed to be the most sensitive to the 
flaring coefficients. 

3- The kinematical behaviour of the warp appeared to be 
dominated by an important velocity gradient. If there is 
indeed a strong warp, so that the HI gas reaches a height 
of a few kpc above the plane at a radius of 20 kpc, then 
we expect a significant rotational velocity gradient per- 
pendicular to the plane. This gradient is of much larger 
amplitude that the simple gradient expected from non- 



cylindrical rotation (appendix A). It is due to the fact 
that the gas orbits in inclined planes (as in the tilted ring 
model of warps) while the bulk of the matter, represented 
by the stars, is flattened in the plane. The periodic orbits 
in the tilted plane are not circular, but elongated, with 
their minor axis in the plane of the galaxy. Their tangen- 
tial velocity is then larger on the minor axis (i.e. in the 
plane) than on their major axis (at their apocenter, at 
their maximum height above the plane). This produces 
a z- velocity gradient of the measured rotational velocity. 
To compute the order of magnitude of this gradient, we 
have taken our model of the NGC 891 potential, including 
the bulge, nuclear disk, and mainly the exponential stellar 
disk. We also introduced the dark halo, with variable flat- 
tening, but this had only a small influence on the orbits 
between 10 and 20 kpc. To find the closed elongated orbits 
in the tilted planes, we then used the shooting method, as 
in the computations of the kinematics of polar rings (e.g. 
Sackett et al 1994, Combes & Arnaboldi 1996). We dis- 
covered that, for all orbits between 10 and 20 kpc average 
radii, the average ellipticity of the orbits is ~ 12%, and the 
tangential velocity on the major axis is ~ 30% less than 
that on the minor axis, for a tilt of the plane of 20° (see 
figure 10). By comparison, the expected non-cylindrical 
z- velocity gradient for the same tilt is of the order of 10% 
only. Assuming such a warp for NGC 891, with a line of 
nodes coinciding with the line of sight, permits us to fit the 
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Fig. 10. Computation of the orbits with the shooting method; 
bottom: shape of closed orbits launched at R=14 and 20 kpc, in 
a plane tilted by 20° from the equatorial plane, the major axis 
is along Ox; top: Time evolution of the radius and velocity of 
the two orbits, one launched at 20 kpc (full line), the other at 
14 kpc (dashed line). The largest relative variations are those 
of the velocity (about 30%) 



observations quite nicely. The HI above the plane on the 
minor axis is taken into account, with a rotational velocity 
less than in the plane, which removes high-Z HI from the 
extreme velocity channels. The high-Z HI is then confined 
in channels closer to the systemic velocity, as shown in 
Fig 8. We can note a posteriori that this effect can produce 
an over-estimation of the rotational velocity at large radii, 
since the velocity in the plane is higher than it should be 
for circular orbits. It is then possible that the true rotation 
curve of NGC 891 is not flat, but slightly falling down at 
large radii. 



The results of this spectral cube modeling are listed below: 

Q W (NGC 891) S3 90° 
H (NGC 891) S3 0.1 (kpc) 

Ci(NGC 891) « 0.02 
Rf(NGC 891) w 12 (kpc) 

C, 2 (NGC 891) R3 0.02 
a Un (NGC 891) S3 0.6 - 0.8 
R W (NGC 891) = 12 (kpc) 

For comparison, we collected the HI parameters for the 
Milky Way, most of them coming from Merrifield (1992): 

B W (MW) S3 80° 
H (MW) = 0.1 (kpc) 

h(MW) S3 0.04 
Rf(MW) S3 16 (kpc) 

( 2 (MW) S3 0.08 
a lin (MW) S3 0.3 
R W (MW) = 12 (kpc) 

We can note significant differences between the two galax- 
ies. Although they both have the same constant HI layer 
thickness and both have a warp beginning at the same 
radius (12 kpc). The warp of NGC 891 is however much 
higher in amplitude (reaching 6 to 8 kpc height at R=22 
kpc, as the Milky Way's warp raises up to 4 kpc height at 
R=25 kpc). Also the flaring of NGC 891 is less by a fac- 
tor 2 in comparison with the Milky Way. We have plotted 
(Fig 11) a perspective view of our best fit model of NGC 
891. 

5. The flattening of an ellipsoidal DMH using the 
outer flaring of the HI layer 

The fit of the rotation curve was made with an ellipsoidal 
dark matter density model. This axisymmetric dark mat- 
ter density model has the advantage to generate an ax- 
isymmetric potential that simplifies substantially the set 
of hydrodynamical and gravitationnal equations formed 
by the hydrostatic equilibrium of the HI gas and the Pois- 
son equation. 

Following Oiling (1995), the projection on the Z axis of 
the Jeans Equation becomes the equation of Hydrostatic 
Equilibrium: 

d(p(Z)al(Z)) 1 d(Rpo\ z ) 

- P(Z)F Z (Z) - dR (9) 

where the term in 9 has been eliminated due to the az- 
imuthal symmetry of the gaseous distribution (m = is 
overwhelming). For Z/R small (Z > 0) the tilt of the 
gaseous velocity dispersion ellipsoid orz is approached by 
(BT87): 

2 / 2 2 \^ 

Vrz ~ (vr - a z)s 
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Fig. 12. We show here the densities of each component of the 
galaxy, for different heights, and for radii corresponding to the 
outer parts. The parameters of these components are found 
from the rotation curve fit. The dotted curve points out the 
HI density vertically distributed (Ho = 0.1 kpc, £ = 0.02) as 
a gaussian following the warp wave (au„(R w = Ylkpc) — 0.8). 
The disk and DMH densities are dashed and full curves re- 
spectively. The DMH has fixed parameters: R c — 14 kpc and 
q — 0.5. The bulge density for R > 10 kpc is totally negligi- 
ble (and do not even appear in the graphs). The dark matter 
density term dominates after 20 kpc. Densities are in units of 
0.001 M Q /pc 3 . 



that we replace in the equation. Although oz cannot 
be measured at large R, observations of face-on external 
galaxies show little or no variation of the Z-dispersion with 
R (e.g. Dickey et al, 1990). Assuming the vertical and ra- 
dial velocity dispersion of the HI gas independent of R: 



2 d\np HI 2 
a? — h (<T R 



2 ^ Z d In phi 



dZ 



galaxy 



R dR 



dz 



(10) 



where 



^galaxy 



bulge 



+ + <P H I + &DMH 



Parallel to the hydrostatic equation (10), let us write the 
Poisson equation for R > 8 kpc (flat region of the rotation 
curve). Deriving with respect to R the equation of the 



Fig. 13. We presents here plots of the right term of the equa- 
tion 16 for different HI parameters (dotted curves). The full 
lines are the left term of eq. 16, DMH densities pdmh, corre- 
sponding to values of q from 0.1 to 0.9. We see that the HI 
term has very different behaviours when we vary the Ho and 
£ parameters, so that the equation can determine "q" unam- 
biguously. Densities are in units of 0.0002 Mq/pc 3 . 



flatness of the rotation curve yields the condition: 



d 2 $galaxy , 1 d$ galaxy 



dR 2 



+ 



R dR 







(11) 



We insert (11) in the Poisson equation written in cylin- 
drical coordinates, which simplifies into (assuming the ax- 
isymmetry of the total galaxy potential): 



AirGp 



galaxy 



galaxy 



dZ 2 



(12) 



We show in Fig 12 the different densities in the outer 
parts, in order to simplify the term 

Pgalaxy = P* + Pbulge + PHI + PDMH 

We see that the bulge density is clearly negligible, while 
after R > 20 kpc, the dark matter dominates largely the 
HI term. We note that the dark matter density considered 
for the graphic corresponds to q = 0.5; it is comparable 
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Fig. 14. We search here the solution of the equation (16) for 
the flaring value of NGC 891: Ci = Ca = C = 0-02. We took here 
ohi isotropic. The value q — 0.2 is a good solution in average, 
the curve being nicely "parallel" to the DMH density curves. 
Densities are in units of 0.0002 M Q /pc 3 . 



for q — 0.8, but twice higher for q = 0.2, so that the 
domination of the dark matter in the outer parts is quite 
certain. We may thus write 



Pgalaxy(R > 20 kpc) ft! PDMH 



(13) 



Combining finally the vertical Jeans equation (10), the 
Poisson equation (12) and the density approximation (13) 
provides: 



Pdmh 



a 2 z d 2 \np HI t r -1 ^ (a R 



2 aV 



4ttG dZ 2 l 4ttG 
d 2 \np HI din phi 



R 



(Z 



dRdZ 



+ 



OR 



)] 



(14) 



The term "In pni(R, Z)" is known since its radial distri- 
bution has been derived in section 2.1 and its vertical dis- 
tribution Zq(R) has been modelled in previous sections: 



Phi(R,@,Z) w pom i exp(- 



(R~R 9 ) 2 
V 



exp(- 



(Z — Z warp (R, 0)) z 
ZZtiR) 



) 



(15) 



The domain of validity of the equation (14) covers R <~ 
20—30 kpc, where the dark matter density dominates. The 



Fig. 15. Solution of the equation (16) for a greater value the 
flaring £ = 0.03, which could have been found without sub- 
stracting the warp. We observe that this value modifies sub- 
stantially the flattening of the dark matter halo, pointing to 
an average value of q ^ 0.4 — 0.5. The method is thus very sen- 
sitive and the allowed range for £ is between 0.015 and 0.035. 
Densities are in units of 0.0002 M & /pc 3 . 



term in bracketts in the right hand side of (14), coming 
from the tilt of the velocity dispersion tensor, showed dur- 
ing computation to be totally negligible (as expected from 
the observations, see Malhotra 1995). For <jr — az = 5 
km/s, this term represents no more than 1 % of the first 
term. So that inserting (15) in the simplified (14) provides 
finally: 



pDMH{pO,DMH{q), R > 20 kpc) = 



4irGZ 2 (R) 



(16) 



We notice that the dark matter density term does not here 
depend on Z, since the information on the flattening "q" is 
in fact included in the central density Po,dmh(q), derived 
from the rotation curve fit. In the equation (16), the verti- 
cal velocity dispersion of the gas, az, carries much of the 
uncertainty. This quantity is frequently approached, in the 
literature, through az ~ <jr, but this is douftful. The in- 
terstellar gas is under the form of small dense clumps, and 
the vertical component of its velocity dispersion should be 
smaller than the radial component, in the same qualita- 
tive manner as the stars. 

We draw in Fig 13 the behaviour of equation (16) in its 
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Fig. 16. Solution of the equation (16) for an anisotropic model 
of the HI velocity dispersion: ohi,z = 7 km/s and o~hi,r, = 10 
km/s. We observe that it rounders substantially the DMH, 
pointing to an average flattening of q ?a 0.4 — 0.5. Densities are 
in units of 0.0002 M Q /pc 3 . 

domain of validity, intending to show its usefulness to de- 
rive q. The figure Fig 14 is the effective solution for our 
problem (i.e. the value of £ = 0.02 found in the previous 
spectral cube modeling section). We see that for a verti- 
cal velocity dispersion equal to the velocity dispersion in 
the plane, the average value of the flattening is q « 0.2, 
pointing to a very flattened DMH. We draw in Fig 15 the 
effects of changing the flaring up to a value of £ = 0.03: 
this rounders the DMH, with q = 0.4 — 0.5. We draw in 
Fig 16 the effects of changing the vertical velocity disper- 
sion of the gas to a lower value of cthi.z = 0.7 <thi,r- This 
makes also the DMH rounder at q = 0.4 — 0.5. 
The method shows up tp be conclusive, reflecting small 
variations of the ^-flaring coefficient and the vertical ve- 
locity dispersion. We also note that there should exist a 
lower limit for the HI vertical velocity dispersion (q < 1.0): 
we find for NGC 891 that <jhi,z > 5.5 km/s. A summary 
of the values of q as a function of flaring parameter and 
HI z- velocity dispersion is listed in Table 2. 

6. Conclusion 

It was unfortunate that NGC 891 was found to have a 
warp which is maximum just along the line of sight, so 



Table 2. Flattening of the DMH of NGC 891 (R c = 14 kpc) 



flaring 


vertical dispersion 


flattening 


c 


ohi,z (km/s) 


q 


0.02 


10 


0.2 


0.03 


10 


0.4-0.5 


0.02 


7 


0.4-0.5 



that the modeling of its HI flaring could not be as ac- 
curate as we could hope. Despite this difficulty, we have 
shown that it is possible to determine the HI plane thick- 
ness and the flaring for this quasi edge-on spiral galaxy 
using sequences of synthetic data cubes. The warp of high 
amplitude for NGC 891 is also a nice explanation for the 
observed strong variation with Z of the rotational veloc- 
ities, a behaviour detected previously by several authors. 
We finally stress that the method seems efficient enough 
to constrain the DMH flattening of every inclined galaxy 
for which we can determine the following quantities: 
&hi,z [from a moment 2 map or face-on galaxies] 
Rc, Po,DMH(q) [from fits of the rotation curve] 
Hq, C [from a 3D spectral cube modeling] 
Noteworthingly, the halo flattening found in this work is 
similar to the result of R. Oiling on NGC 4244, both sug- 
gesting very flattened halos. It thus may become now in- 
teresting to test the dark matter model predicted by Pfen- 
nigcr & Combes, a warped disk of H2 cold molecular gas. 
We shall present in a forthcoming paper the results of our 
method applied to other galaxies with interferometric HI 
data cubes. 
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7. Appendix A: Non-cylindrical motions in galax- 
ies 

Let us consider a very flattened Kuzmin potential to rep- 
resent the whole galaxy, 



*{R,Z) = 



GM 



^R 2 + (a + \Z\) 2 



We expand v 2 — R^ assuming a <C Z <C R, obtaining: 

Vc(z) = «0(1-^(|) 2 ) 

Applying this formula for NGC 891, at R = 10 kpc and 
Z = 3 kpc gives for u = 225 km/s a loss of velocity 
of about 15 km/s. A vertical gradient of velocity is thus 
generally expected in galaxies. We include explicitly this 
gradient in the model in multiplying the disk velocity by a 
factor fv(Z), which is chosen in practice to be (1^|(^) 2 )- 

8. Appendix B: Derivation of the spiral arm veloc- 
ities 

We entered in the modeling of the spiral arms their per- 
turbed velocities associated with the mode to (practically, 
m = 2). Using Euler's equation, BT87 estimated that Vr 
and Vq are of the same order and thus the continuity equa- 
tion could be simplified, (pages 356-359) resulting in the 



formula of the perturbed radial velocity associated with 
the mode to: (under adaptations for the present context) 

V R (m) = -?*h^L ^ [I C os(A :mJ R + m(e-e m ))] 

The perturbed tangential velocity associated with the 
mode to is obstained from formula (6-36) of BT87, yield- 
ing: 



Ve(m) = 



-IB A r< 



1 



sin(fc m i? + to(9 - m ))] 



k m A R 

"B" being Oort's second "constant". It is noteworthy 
that Vr changes sign at corotation, while V@ does not. 
The kincmatical parameters a-priori unknown in the for- 
mula above are derivable from a multi-component fit of 
the rotation curve (Fig 4) : the value of = -^r 1 - , k = 



R 



dR 



4Q, 2 and B = B(R) = -fw. The main modifi- 



cation in the spectral cube construction due to the con- 
sideration of these perturbed velocities, is to change the 
formula of the line-of-sight velocity which is now: 

V los (m = 0) = V sys + V rot sin I cos Qf v (Z) 
Vi os (m^0) = V sys + [(Vrot + Ve(m)) sin I cos 6 
+V R (m)smIsmQ]f v (Z) 

9. Appendix C: potentials for generalized dark 
matter densities 

Simulations of spherical infall in an Einstein-de Sitter uni- 
verse onto seeds perturbations and onto density peaks 
in power cosmologies result in objects with p cx Rr a , 
1.6 < a < 2.25 (Dubinski & Carlberg, 1991). This is close 
enough to the isothermal profile (i.e. R~ 2 ) 7 which is of- 
ten adopted for the sake of simplicity. However, Lake & 
Feinswog (1989) showed that profiles in R~ 3 and even i?~ 4 
were also possible to obtain a good fit of observed rota- 
tion curves. Taking this result into account, we consider 
hereafter the generalized family of densities: 
Po,dmh 



(P>0) 



This density model recovers, for p = 1, the formula of the 
isothermal density profile used in this article and previ- 
ously by other authors. We give here the formula of the 
dark matter potential for p =/= 1, which can be derived us- 
ing formula 2.6 of de Zeeuw & Pfenniger (1988), and the 



new variable (x 



DMH 



(R,Z) 



^Ju+qlRl 
27T 



), providing: 



G Po 



DMH 



qR 2 c 



I 



1-p 



dx 



,0 e 2 x 2 + 1 

The rotation curve contribution is, from that moment, 
straightforward. 
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